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Modeling  Fusion  of  Cellular  Aggregates  in  Biofabrication  Using 

Phase  Field  Theories 

Xiaofeng  Yang*  Vladimir  Mironov'and  Qi  Wang* 


Abstract 

A  mathematical  model  based  on  a  phase  field  formulation  is  developed  to  study  fusion  of 
cellular  aggregates/clusters.  In  a  novel  biofabrication  process  known  as  bioprinting  [25],  live  mul¬ 
ticellular  aggregates/clusters  are  used  to  make  tissue  or  organ  constructs  via  the  layer-by-layer 
deposition  technique  in  compatible  hydrogels  rich  in  maturogen;  the  bio-constructs  embedded 
in  hydrogels  are  then  placed  in  bioreactors  to  undergo  the  fusion  process  of  self-assembly,  matu¬ 
ration,  and  differentiation  to  form  the  desired  functional  tissue  or  organ  products.  We  formulate 
the  mathematical  model  to  study  the  morphological  development  of  the  printed  bio-constructs 
during  fusion  by  exploring  the  chemical-mechanical  interaction  between  cellular  aggregates  in¬ 
volved.  Specifically,  we  treat  the  cellular  aggregates  and  the  surrounding  hydrogels  as  two 
immiscible  complex  fluids  and  then  develop  an  effective  mean-field  potential  that  incorporates 
the  long-range,  attractive  interaction  between  cells  as  well  as  the  short-range,  repulsive  inter¬ 
action  due  to  immiscibility  between  the  cell  and  the  hydrogel.  We  then  implement  the  model 
using  a  high  order  spectral  method  to  simulate  the  making  of  a  set  of  tissues/organs  in  sim¬ 
ple  geometries  like  a  ring  or  a  sheet  of  tissues  and  a  Y-  or  T-shaped  vascular  junction  by  the 
layer-by-layer  deposition  of  spheroidal  cellular  clusters  in  the  bioprinting  technology. 


1  Introduction 

Tissue  fusion  is  an  ubiquitous  phenomenon  during  embryonic  development  and  morphogenesis 
[30].  The  natural  tissue  fusion  in  vivo  usually  occurs  in  two  steps.  It  starts  from  the  initial  tissue 
opposition  and  follows  by  sequential  actual  tissue  fusion  after  establishing  direct  contacts  between 
adjacent  embryonic  tissues.  The  impaired  tissue  fusion  process  during  embryonic  development 
resulted  in  embryonic  malformation  and  defects  for  example  such  as  cleft  palate.  In  the  past, 
bioengineering  processes  have  been  devised  to  fabricate  tissues  under  controlled  conditions  using 
tissue  self-assembly,  in  which  one  seeds  cells  into  biodegradable  polymer  scaffolds  or  gels,  which 
are  then  cultured  in  bioreactors  for  several  weeks  and  finally  implanted  into  the  recipient  organism, 
where  the  maturation  of  the  new  organ  takes  place  [17,  24], 

In  a  novel  biomimetic  biofabrication  process,  called  “bioprinting” ,  multicellular  tissue  spheroids 
or  aggregates  are  used  as  fundamental  building  blocks  to  construct  the  3-D  tissue  or  organ  [12,  17, 
24,  25].  The  multicellular  aggregates  are  first  prepared  in  the  form  of  tissue  spheroids  that  consist 
of  thousands  of  cells  blended  with  bio-compatible  hydrogels.  They  are  then  directly  deposited 
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by  computer-aided  design  tools  into  desired  3-D  tissue  or  organ  constructs  via  the  layer-by-layer 
deposition  technique.  The  bio-constructs  immersed  in  a  hydrogel  are  then  placed  in  bio-incubator 
for  maturation.  During  the  incubating  process,  the  cellular  aggregates  immersed  in  compatible 
hydrogel  are  expected  to  fuse  into  a  3-D  tissue  or  organ  following  the  natural  rule  of  histogenesis 
and  organogenesis  [17,  29,  18,  19].  During  this  process,  tissue  self-assembly,  fusion,  differentiation 
and  maturation  can  occur.  Experimental  evidence  demonstrates  that  setting  lumenized  vascular 
tissue  spheroids  in  a  roll,  as  the  result  of  tissue  fusion,  a  lumenized  vascular  tube  can  be  formed  [25]. 
Tissue  fusion  driven  is  a  fundamental  biophysical  process  in  emerging  organ  bioprinting  technology. 

The  bio-constructs  ranging  from  the  ones  comprised  of  tissue  spheroids  to  functioning  tissues 
or  organs  all  exhibit  fluid  behavior  during  tissue  fusion  processes  and  are  bona  fide  soft  materials 
with  various  degree  of  viscoelasticity.  Numerical  simulations  of  the  morphogenesis  phenomenon 
in  bioprinting  using  the  Monte  Carlo  method  and  experimental  evidences  [17,  18,  19,  29,  28,  24] 
pointed  out  the  strong  influence  of  surface  tension  to  the  tissue  spheroids  fusion  in  the  biofabrication 
process  when  the  biomaterials  (tissue  spheroids  and  hydrogels)  are  regarded  as  viscous  fluids.  A 
separate  study  using  experimental  methods  and  numerical  simulation  with  the  cellular  Potts  model 
also  explored  the  effect  of  surface  tension  and  the  hydrodynamics  to  the  rounding  of  cell  aggregate 
[27].  The  experimental  evidence  and  the  Monte  Carlo  simulation  also  demonstrate  an  interesting 
compaction  phenomenon  observed  initially  in  the  self-assembly  of  the  bio-construct  before  it  evolves 
into  a  smooth  tissue  or  organ  [24,  26,  41,  32],  This  phenomenon  resembles  the  long  rang  attractive 
interaction  among  small  particles,  but  it  occurs  at  a  much  large  length  and  time  scale.  The 
molecular  origin  of  the  compaction  can  be  the  consequence  of  cell  motility  or  chemically  induced 
electrostatic  interaction  between  cells  mediated  by  the  hydrogel.  At  large  length  and  time  scales, 
the  phenomenon  can  be  modeled  using  a  coarse-grain  model  featuring  the  long  range  aggregate- 
aggregate  interaction,  which  will  be  detailed  in  the  next  section. 

It  is  known  that  cell  aggregates  exhibit  fluid-like  behavior  during  the  tissue  fusion  process 
at  large  length  and  time  scales  [24],  Both  the  multicellular  cell  aggregates  in  the  form  of  tissue 
spheroids  and  the  hydrogel  in  which  the  tissue  spheroids  are  embedded  are  multiphase  complex 
fluids  in  nature.  Together,  they  form  a  binary  complex  fluid  system,  in  which  one  fluid  is  the 
cellular  aggregates  and  the  other  is  the  hydrogel  at  the  coarse-grain  level.  The  dynamics  of  the 
binary  fluid  system  can  then  be  modeled  in  the  coarse-grain  length  scale  by  models  suitable  for 
multiphase  complex  fluids.  Given  the  large  time  scale  in  tissue  fusion,  elasticity  of  the  tissue 
spheroids  and  the  bioconstructs  made  up  of  the  tissue  spheroids  can  be  effectively  ignored  unless 
an  imposed  external  field  is  introduced.  We  therefore,  model  the  bioconstructs  at  this  time  scale 
as  viscous  fluids. 

An  excellent  review  for  biofabrication  methods,  physical  mechanisms  and  the  discussion  on  the 
fluidity  cellular  aggregates  is  available  in  [24] .  An  elasto- visco-plastic  continuum  model  is  developed 
to  investigate  the  cell  aggregate  deformation  properties  in  tension  [31],  in  which  a  trigger  for  yield 
stress  is  implemented.  In  large  time  scale,  the  model  reduces  to  the  viscous  limit  as  well.  The 
biophysical  mechanism  of  tissue  fusion  process  is  not  completely  understood.  Modern  mathematical 
modeling  and  computer  simulation  opens  unique  opportunities  for  explaining  physical  nature  of 
tissue  fusion  process  and  predicting  possible  undesirable  outcomes  which  is  essential  for  designing 
optimal  bioprinting  protocol. 

Modeling  and  simulating  immiscible  multiphase  fluid  flows  have  been  challenging  both  mathe¬ 
matically  and  technically  over  the  years.  Various  mathematical  theories  and  computational  tech¬ 
nologies  have  been  developed  to  tackle  the  problem.  The  front  tracking  method  [10],  boundary 
integral  method  [16],  level-set  method  [35,  33],  volume-of-fluid  method  [15],  immersed-boundary 
method  [34,  23],  and  phase  field  method  [1,  38,  37,  36]  have  all  been  proposed,  implemented  and 
refined,  each  of  which  has  shown  effectiveness  in  designated  applications.  Some  of  the  methods 
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are  comparable  in  designing  concept  while  others  can  be  combined  to  yield  more  effective  com¬ 
putational  technology  [22],  Among  all  above,  the  phase  field  method  for  multiphase  fluid  flows  is 
perhaps  the  simplest  to  implement,  given  that  the  phase  boundary  is  embedded  in  a  level  set  of  the 
phase  variable  governed  by  a  dissipative  evolutionary  equation,  and  the  most  physically  relevant. 
Because  of  the  ease  of  use,  simplicity  and  relevance  to  material’s  physics,  refined  details  in  the 
formulation  of  the  phase  field  equation  can  be  carefully  carved  out  and  the  free  energy,  especially, 
the  interfacial  free  energy  for  the  multiphase  fluid  can  be  devised  properly  to  ensure  accuracy  in 
numerical  computations  and  fidelity  in  physical  modeling. 

In  an  immiscible  binary  fluid,  the  phase  field  method  employs  a  phase  variable  0  <  0  <  1  to 
track  each  phase  in  the  binary  fluid:  0  =  1  describes  the  phase  of  fluid  I  and  0  =  0  denotes  the 
one  occupied  by  fluid  II  while  0  <  0  <  1  describes  the  interfacial  region.  The  phase  variable  is  also 
known  as  a  labeling  function  for  identification  purposes.  The  time  evolution  of  the  phase  variable 
0  is  governed  according  to  the  Cahn-Hilliard  equation  [4,  5] 

'[]’  V  •  1 AV //  j.  (1.1) 

where  A  is  proportional  to  the  mobility  and  n  is  the  chemical  potential  of  the  multiphase  fluid 
system,  a  functional  of  the  phase  variable  0,  and  ^  is  the  material  derivative.  The  phase  variable 
0  can  be  identified  with  the  volume  fraction  of  fluid  I;  so,  1  —  0  serves  as  the  volume  fraction  of  fluid 
II.  In  this  formulation  of  the  transport  equation  for  0,  the  flux  of  0  is  assumed  to  be  proportional 
to  the  force  due  to  the  prescribed  chemical  potential.  However,  a  more  physically  appropriate 
assumption  is  to  assume  the  excessive  transporting  velocity  of  0,  in  addition  to  the  bulk/average 
velocity,  is  proportional  to  the  force  due  to  the  chemical  potential,  a  consequence  of  the  friction 
dynamics  [8,  3].  This  leads  to  the  singular  or  modified  Cahn-Hilliard  equation 

d('n  V  •  (A S0VM),  (1.2) 

where  As  is  the  mobility.  If  we  assume  mixing  only  takes  place  in  the  interfacial  layer,  it  would  be 
reasonable  to  assume  As  =  A],(l  —  0),  where  A),  is  a  constant.  Hence,  the  volume  fraction  will  retain 
a  constant  within  the  bulk  of  each  fluid.  In  the  classical  Cahn-Hilliard  equation,  A  is  a  constant; 
whereas  it  is  a  phase  variable  dependent  function  in  the  modified  case.  As  A  — >•  0,  the  transport 
equation  for  the  phase  variable  reduces  to  a  pure  transport  equation  for  the  phase  boundary  which 
is  commonly  used  in  the  level-set  method  [35,  33]. 

The  phase  field  model,  also  known  as  the  diffuse  interface  method,  benefits  from  the  dissipation 
mechanism  in  the  transport  equation  for  0.  There  exists  an  interfacial  layer  in  which  two  fluids  mix 
due  to  dissipation.  It  is  expected  to  maintain  a  constant  value  for  the  phase  variable  0  outside  the 
layer.  For  a  slow  varying  interface,  namely,  the  interface  with  a  small  or  intermediate  curvature, 
phase  field  models  yield  acceptable  interfaces  with  thickness  controlled  within  a  few  grid  points. 
When  the  curvature  becomes  large,  resolving  the  interfacial  layer  accurately  becomes  a  challenge. 
Nonetheless,  the  material  system  consisting  of  cellular  aggregates  and  the  host  hydrogel  matrix  does 
not  necessarily  maintain  a  sharp  interface  at  the  coarse-grain  level  given  the  cell  motility.  Hence, 
the  phase  field  method  may  serve  as  a  good  systematic  way  to  model  this  biomaterial  system. 

To  formulate  the  cellular  spheroid  fusion  in  a  hydrogel  matrix,  we  identify  the  cellular  spheroid 
as  fluid  I  and  the  ambient  fluid,  the  hydrogel,  as  fluid  II  in  the  framework  of  phase  field  theories. 
We  account  for  three  distinct  interaction  forces  between  these  two  “immiscible”  fluids.  First, 
we  introduce  a  short-range  repulsive  bulk  potential  to  maintain  the  clear  separation  between  the 
two  distinctive  fluids,  a  force  to  maintain  immiscibility  between  the  two  fluids;  second,  we  use  a 
configurational  entropic  term  in  the  interaction  potential  to  describe  the  philic  interaction  between 
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the  cellular  spheroids  in  close  proximity  and  the  resultant  surface  tension  when  coupled  with  the 
bulk  potential;  third,  we  incorporate  a  long-range  attractive  interaction  potential  based  on  collective 
Lennard- Jones  interaction  among  the  cells  physically.  The  third  force  is  necessary  to  capture  the 
compaction  process  observed  in  the  initial  self-assembly  process  of  cellular  aggregates.  Our  phase 
held  model  is  primarily  developed  based  on  these  interaction  forces  and  their  coupling  with  the 
mass  and  momentum  transport  of  the  binary  fluid  system.  Should  additional  forces  between  the 
approaching  cellular  spheroids  be  identified  experimentally  at  the  cellular  level,  a  potential  yielding 
the  mean  force  can  be  derived  accordingly. 

We  organize  the  rest  of  the  papers  into  three  sections.  In  the  second  section,  we  derive  the 
detailed  phase  held  model  and  discuss  its  nondomensionalization.  In  the  third  section,  we  design 
an  efficient  numerical  scheme  to  discretize  the  governing  system  of  partial  differential  equations 
and  discuss  the  numerical  implementation  that  simulates  the  cellular  cluster  fusion  process  to 
form  tissues  and  organs.  Tissues  in  simple  geometries  like  rings,  sheets  are  simulated  along  with 
biofabrication  of  bifurcating  Y  or  T-shaped  vascular  veins  via  layer-by-layer  deposition  in  the  fourth 
section. 


2  Mathematical  formulation  of  the  phase  field  model 


We  treat  the  cellular  cluster  as  a  blob  of  complex  fluids  and  its  surrounding  fluid  (hydrogel 
matrix)  as  a  viscous  fluid.  The  cellular  cluster  along  with  the  surrounding  huid  then  forms  a 
binary  huid  mixture  of  two  immiscible  fluids.  The  binary  huid  system  is  assumed  incompressible, 
which  is  a  reasonable  approximation.  We  use  a  phase  variable  4>  to  label  each  huid  phase  in  the 
system: 


f  1,  in  cellular  cluster, 

1  0,  in  host  huid  matrix. 


(2.1) 


We  denote  the  average  velocity  of  the  huid  mixture  by  u.  The  transport  equation  for  the  mass  and 
momentum  of  the  mixture  system  is  governed  respectively  by 


Pw  =  V  '  (<M  +  (!  -  (t))T2)  -  (Vp  +  4>Vp) 
V  •  u  =  0, , 


(2.2) 


where  p  =  4>pi  +  (l  —  (j))p2  is  the  effective  density  for  the  binary  huid,  pi,  t\  and  p2,  T2  are  the  density 
and  the  extra  stress  tensor  for  the  huid  consisting  of  cellular  clusters  and  the  host  huid  matrix, 
respectively,  p  is  the  extended  chemical  potential  of  the  mixture  system,  and  p  is  the  hydrostatic 
pressure.  The  chemical  potential  of  the  material  system  is  calculated  from  the  system  free  energy 
consisting  of  the  extended  “mixing  free  energy”  [9,  13]  defined  by 

Fmix{(t>)  =  J  f  ^ikT\VcP\2  +  72 kTF(</>))dx,  (2.3) 

z  Jn 

where  fl  is  the  domain  that  the  mixture  occupies,  F(<j))  =  02(1  —  (j))2  is  the  Ginzburg-Laudau  double 
well  potential,  k  is  the  Boltzmann  constant,  T  is  the  temperature,  71  is  a  parameter  measuring  the 
strength  of  the  conformation  entropy  and  72  is  the  strength  of  the  bulk  mixing  free  energy.  We 
introduce  an  interaction  potential  due  to  the  long-range  cellular  interaction: 

FCei{(t>)  =  J ^  j^x{\x  ~  y\)4>{x,t)(j){y,t)dy  \  dx.  (2.4) 
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We  assume  that  the  long-range  cellular  cluster  interaction  is  due  to  the  collective  attractive  inter¬ 
action  between  the  cells  parametrized  by  the  Lennard-Jones  potential  y(r)  with 


X{r)  =  Vlj(t)  =  4e  ((^)12  -  (^)6)  , 


(2.5) 


where  e  is  a  suitably  chosen  energy  constant  which  is  influenced  by  the  cellular  environment  around 
the  surface  of  the  cellular  cluster  (i.e.,  it  can  be  a  feedback  function  should  more  detailed  signal¬ 
ing  molecular  information  becomes  available  when  cellular  clusters  are  placed  in  proximity),  a  is 
the  finite  distance  at  which  the  inter-particle  potential  is  zero  and  r  is  the  distance  between  the 
interacting  “particles”  (cells)  [20].  Like  e,  a  is  a  parameter  that  depends  on  the  detail  of  cellular 
interaction  when  used  in  this  system.  Thus,  the  total  free  energy  of  the  system  is  then  given  by 


Ftot.  —  Fn 


+  Fcei . 


(2.6) 


The  extended  chemical  potential  for  the  material  system  is  calculated  by  the  variational  deriva¬ 


tive 


5Ft 


V  = 


tot 


=  -71 V  (j>  +  72(1  —  </>)</>(l  —  20)  +  /  x{\x  ~  y\)(j)(y,t)dy. 
$(P  Jn 


In  practice,  we  use  the  truncated  Lennard-Jones  potential 


X(r)  = 


VLj{r)  -  VLJ{rc)  r  <  rc, 
0  r  >  rc, 


(2.7) 


(2.8) 


where  rc  >  2.5 a.  In  the  calculations  we  conducted  in  this  paper,  we  use  rc  =  2.5c.  We  remark  that 
the  Lennard-Jones  potential  is  attractive  at  the  “long-range”  and  repulsive  in  the  “short-range”; 
its  use  in  this  model  is  primarily  for  the  “long-range”  attractive  effect.  The  short-range  repulsive 
effect  of  the  Lennard-Jones  potential  is  easily  offset  by  the  conformational  entropy  which  promotes 
a  different  kind  philic  interaction  at  the  short  to  mid-range  interaction. 

The  phase  variable  4>  or  the  volume  fraction  of  the  cellular  component  is  transported  via  the 
Cahn-Hilliard  equation 


^  +  V-(u^)  =  V-A(V/x), 


(2.9) 


where  A  is  the  mobility  that  is  normally  a  function  of  the  volume  fraction  cj).  Often,  it  is  replaced 
by  a  constant  for  simplicity.  At  the  interface  between  the  two  fluid  phases,  from  the  Cahn-Hilliard 
equation,  we  identify  the  velocity  of  the  cellular  fluid  as 


ui  =  u  —  —  V/x. 


(2.10) 


and  that  of  the  host  fluid  matrix  as 


u2  =  u  + 


1  - 


:  V/J. 


(2.11) 


These  two  velocities  are  identifiable  only  within  the  interfacial  layer  between  distinctive  phases. 
Within  each  phase,  we  assume  the  complex  fluid  homogeneous. 

The  extra  stress  tensors  are  given  according  to  the  material’s  property  of  each  fluid  phase  in¬ 
volved.  Even  though  the  multicelluar  cluster  are  Non-Newtonian  in  nature,  the  rheological  response 


5 


in  long  time  is  approximately  viscous  or  Newtonian.  Considering  the  large  time  scale  involved  in 
the  morphological  development  process,  we  assume  both  fluids  viscous  with  distinct  viscosities: 


n  =  2771 D,  r2  =  27/2D, 


(2.12) 


where  771 , 772  are  the  viscosity  for  the  cellular  cluster  fluid  and  the  host  fluid  matrix,  respectively,  D 
is  the  rate  of  strain  tensor  for  the  mixture  associated  to  the  average  velocity  field  v  defined  by 

D  =  1(Vu  +  (Vu)t). 

We  investigate  the  interfacial  dynamics  of  the  binary  fluid  associated  with  the  cellular  cluster 
fusion  in  a  2-D  domain  fl  =  [0,  H]2  and  3-D  domain  fl  =  [0,  H]3.  At  the  boundary  of  the  computed 
domain  <90,  we  impose  no- flux  boundary  conditions  for  the  phase  variable  of  the  cellular  fluid,  the 
Dirichlet  boundary  condition  for  the  velocity: 

V0  •  n|an  =  0, 

<  (u</>  -  AVp)  •  n|an  =  0,  (2.13) 

,  u|an  =  0. 


We  use  a  characteristic  time  scale  to  and  a  length  scale  h  to  nondimensionalize  the  variables 

7  t  _  X  _  Uf0  _  pt% 

t  =  T0’x=h’u=lT’p  =  ^-  <2'14) 


The  length  scale  h  is  determined  by  the  computational  geometry  while  the  time  scale  is  done  by 
either  the  growth  time  scale  of  the  interface  or  the  cellular  cluster  fusion  time  scale.  The  following 
dimensionless  equations  then  arise 


A 


'yikTt-Q  -p 

PohA  1  12 


l2kTtl 
Poh 2  ’ 


Re2 


_  Poh2 
V2to  ’ 


Re  1  = 


Poh2 
Vito  ’ 


P  =  ^  +  (1_0) 


poh2  ’ 


(2.15) 


where  Re ±  and  i?e2  are  the  Reynolds  number  for  the  cellular  fluid  and  the  ambient  fluid,  respectively, 
po  is  an  average  (or  a  constant  reference)  density.  For  simplicity,  we  drop  the ~ on  the  dimensionless 
variables  and  the  parameters.  The  system  of  governing  equations  for  the  binary  fluid  in  these 
dimensionless  variables  are  given  by 


p—  =  V  •  {<j>Ti  +  (1  -  4>)t2)  -  (Vp  +  </>Vp) 

V  •  u  =  0 

^  +  V-(<H  =  V-(AVp) 

where  t\  =  D,  r2  =  The  dimensionless  energy  density  is  now  given  by 

/  =  yl  Wl2  +  T2<Kl  -  </>)(!  -  20)  +  ^  x(\x  -  y\)(j)(y,  t)dy. 


(2.16) 

(2.17) 

(2.18) 


(2.19) 
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3  Numerical  Methods 


We  apply  the  phase  field  model  developed  in  the  previous  section  to  study  the  fusion  of  the 
cellular  clusters  when  they  are  embedded  in  a  host  hydrogel  matrix,  i.e.,  the  time  evolution  of 
the  overall  exterior  boundary  of  the  cellular  clusters  arranged  in  a  designer’s  pattern  for  tissue 
and  organ  generation.  We  first  describe  the  numerical  scheme  to  solve  the  Cahn-Hilliad  equation 
without  the  velocity  (2.18),  i.e.  zero  velocity  field.  This  decoupled  system  is  used  to  simulate  the 
fusion  of  spheroids  which  are  in  contact  initially.  For  this  case,  the  bulk  velocity  remains  near  zero 
when  fusion  starts.  Thus  we  only  need  to  focus  on  the  binary  fluid  evolution  without  significant 
fluid  motion  and  employ  the  Cahn-Hilliard  equation  to  investigate  the  interfacial  dynamics  at  this 
moment.  Furthermore,  we  will  discuss  the  numerical  scheme  to  solve  the  coupled  system  of  Cahn- 
Hilliad  and  Naviei'-Stokes  equations  (2.16-2.18)  to  investigate  the  hydrodynamic  effect  that  drives 
cellular  spheroid  fusion  when  they  away  from  each  other  initially. 

We  denote  T  =  AFi  and  ? j  =  The  Cahn-Hilliard  equation  takes  the  following  form: 


+  rA(A<£  -  (/'(</>)  +  fi))  =  0, 


(3.1) 


where  f{$)  =  4^(1  -  4>){l  -  2<£)<£,  ft  =  ±  fQ  x(\x  -  y\)<f>(y,t)dy. 

We  use  the  first-order  backward  Euler  method  to  discretize  the  time  derivative.  The  numerical 
scheme  reads: 


vn+l 


TS 


St 


+ r  AAr+1  -  ~^rr+l  -  r)  =  pa  (f\r)  +  Mr)), 


(3.2) 


where  the  bilaplace  term  is  treated  implicitly  and  all  other  nonlinear  terms  are  treated  explicitly. 
An  extra  term  associated  with  the  artificial  parameter  S  is  introduced  in  order  to  balance  the 
unstable  constraint  condition  on  the  time  step  due  to  the  explicit  treatment  of  nonlinear  terms 
[39].  The  spatial  discretization  of  the  semi-discrete  equation  is  carried  out  using  spectral-Galerkin 
method  with  the  Legendre-Gauss-Lobatoo  points  in  one  direction  and  Fourier  method  in  the  other 
two  directions.  Hence,  in  our  implementations,  the  boundary  condition  in  one  direction  is  physical 
and  the  other  two  directions  are  periodic. 

For  the  coupled  system  of  Cahn-Hilliard  equation  and  the  Navier-Stokes  equations  (2.16-2.18), 
we  use  the  first-order  pressure-correction  scheme  developed  in  [6,  7,  40,  11].  For  simplicity,  we 
assume  the  density  of  the  cellular  fluid  and  the  hydrogel  is  the  same  and  then  set  the  density  p  =  1 
in  our  nondimensionalization.  Assuming  that  ( un ,  un,pn ,  <f>n)  is  known,  the  numerical  scheme  reads: 


1.  Update  of  the  intermediate  velocity  field  un+1: 

1  1 


un+i  -  un 


St 


-  (fHW  +  FHw)A(un+1  -  u")  +  Vpn  =  N(r,  u"), 


'  2Re\  2 Re2 


(3.3) 


un+1|an  =  0, 


where  N(</>,  u)  =  -u  •  Vu  -  {fVp  +  (7^  +  kf  “  2ifel  “  2iUv2u- 
2.  Enforcing  the  divergence-free  condition  through  projection  to  obtain  velocity  and  pressure 

(■ un+\pn+ 1); 


u”+l  _  ■Qn+1 


+  V  (pn+1  —  pn)  =  0, 


St 

d(pn+ 1  —  pn) ! 
dn 


Ian 


=  0, 


(3.4) 
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3.  Update  the  phase  variable  (j)n+1: 


v  +  un+1  •  vr  +  TAA<pn+i  -  -4 A(r+1  -  <n  =  rA u'm  +  mo, 

Ot  T]z 

d^n+l  dA  cj)n+1 

— o — an  =  0)  — x - an  =  0, 

on  on 


where  n  is  the  outward  normal  on  the  boundary.  We  implemented  the  coupled  system  in  2-D  and 
the  decoupled  system  in  full  3-D.  In  the  2-D  numerical  results  presented  below,  we  use  129  Legend 
basis  in  one  direction  and  129  Fourier  modes  in  the  other  direction.  In  the  3-D  simulations,  we  use 
129  Fourier  modes  for  the  extra  direction.  The  computed  domain  is  [—1, 1]  x  [—1, 1]  for  the  2-D 
case  and  [—1, 1]  x  [—1, 1]  x  [—1, 1]  for  the  3-D  case.  For  the  2-D  coupled  model  (2.16-2.18),  we  set 
the  parameters  in  the  momentum  equations  as  Re \  =  Re 2  =  10-10,A  =  10  .  In  all  simulations, 
we  always  fix  T  =  0.0002, 77  =  0.01,  S  =  2,  rc  =  0.5,  e  =  0.0005,  a  =  0.2. 


4  Results  and  Discussions 

We  next  investigate  cellular  spheroid  fusion  with  the  numerical  code  developed  in  the  previous 
section.  We  first  look  into  the  compaction  phenomenon  in  bioprinting  of  tissues  using  the  2-D  Cahn- 
Hilliard/Navier-Stokes  solver.  In  the  bioprinting  process,  when  the  cellular  spheroids  are  deposited 
into  a  designed  tissue  construct,  the  cellular  clusters  undergo  an  initial  compaction  before  they  start 
fusing.  During  the  compaction  process,  the  center  of  mass  of  each  cellular  spheroid  moves  towards 
the  geocenter  of  the  construct  under  an  attractive  force.  This  motion  occurs  on  a  longer  time  scale 
than  the  fusion  process.  As  an  illustration,  we  will  simulate  the  evolution  of  a  ring  shaped  bio¬ 
construct  during  the  biofabrication  process.  Then,  we  simulate  the  fusion  process  for  the  cellular 
spheroids  already  in  contact  using  the  3-D  Cahn-Hilliard  solver  for  a  set  of  bio-constructs  in  simple 
geometries. 

4.1  Initial  attractive  motion  of  the  cellular  spheroids  in  the  tissue  ring  forma¬ 
tion 

The  bio-construct  made  of  the  cellular  spheroids  with  the  spheroids  initially  positioned  slightly 
apart  undergoes  a  structural  contraction  globally.  This  is  shown  experimentally  in  [17]  as  the 
compaction  phenomenon.  We  first  examine  this  phenomenon  using  the  2-D  Navier-Stokes/Cahn- 
Hilliard  solver.  We  investigate  how  the  attractive  motion  of  the  cellular  spheroids  deposited  in 
a  ring  pattern  evolves  in  time.  The  lay  out  of  the  ring  construct  consists  of  8  cellular  spheroids 
positioned  slightly  apart.  Figure  1  depicts  the  configurational  layout  of  the  8  cellular  spheroids  at 
different  time  of  the  simulation.  At  the  initial  stage  of  fusion,  the  center  of  mass  of  each  spheroid 
moves  towards  the  geo-center  of  the  system  under  the  long-range  attractive  interaction  among  the 
cells  leading  to  an  evolved  ring  pattern  with  the  cellular  spheroids  in  contact  to  each  other.  As  soon 
as  the  spheroids  get  in  contact  with  one  another  the  philic  interaction  due  to  the  conformational 
entropy  takes  over.  A  smooth  tissue  ring  consisting  of  the  cells  forms  under  the  dominating  impact 
of  the  hydrodynamic  surface  tension.  Given  the  fact  that  the  volume  of  the  cellular  region  is 
conserved  in  the  model,  the  radius  of  the  smoothed  ring  keeps  shrinking.  The  radius  reduction  in 
time  until  the  spheroids  are  in  contact  is  recorded  in  Figure  2,  which  agrees  qualitatively  with  [26]. 
The  continuing  ring  shrinking  phenomenon  during  fusion  after  the  cellular  spheroids  are  in  contact 
will  be  discussed  in  more  details  next. 


Figure  1:  Fusion  of  a  ring  construct  from  ten  deposited  spheroidal  cellular  clusters.  Initially, 
the  ten  spheroids  are  placed  equidistance  to  the  geocenter  and  completely  separated  from  each 
other.  The  long  range  attractive  force  moves  the  center  of  mass  of  each  spheroid  towards  the 
geocenter.  Later,  surface  tension  takes  over  to  facilitate  the  fusion  process.  The  snapshots  are  for 
t  =  5, 40, 55, 57, 58, 60, 80, 85. 


Figure  2:  The  external  radius  of  the  ring  construct  as  a  function  of  time  up  to  t  =  60.  Compaction 
of  the  ring  construct  and  later  contraction  of  the  smooth  ring  is  evident. 
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Figure  3:  Fusion  of  a  ring  construct  from  eight  deposited  spheroidal  cellular  clusters.  Surface 
tension  plays  a  major  role  in  the  fusion  process.  The  snapshots  are  taken  at  t  =  0.1, 2, 16, 30, 38. 


This  simulation  demonstrates  the  effect  of  the  long-range  attractive  interaction  among  the  cells 
to  the  initial  configurational  evolution  of  the  cellular  spheroids  within  a  designed  bio-construct.  For 
closely  packed  spheroids  that  are  already  in  contact  this  effect  is  secondary.  The  primary  force  in 
that  situation  is  the  philic  conformational  entropic  one  closely  related  to  the  surface  tension  effect. 
Therefore,  in  the  following  simulation,  we  will  instead  focus  on  the  interaction  dominated  by  the 
surface  tension  in  the  spheroid  fusion  process  upon  contact  in  3-D. 

4.2  Morphological  evolution  of  a  ring  construct  made  of  cellular  spheroids 

Using  the  3-D  Cahn-Hilliard  solver,  we  investigate  how  a  tissue  ring  can  be  fabricated  by 
depositing  a  ring  of  cellular  spheroids  in  contact  in  the  bio-printing  technology  and  then  watch 
how  this  ring  can  eventually  contract  into  a  sphere  over  a  long  period  of  fusion  process.  The  initial 
deposited  cellular  spheroids  in  a  construct  of  the  ring  shape  consists  of  eight  spheroidal  cell  clusters. 
Due  primarily  to  the  philic  effect  (dissipation  effect)  initially,  the  spheroids  join  into  a  corrugated 
ring;  then,  the  surface  tension  takes  over,  the  surface  is  gradually  smoothed  out.  In  the  meantime 
the  radius  of  the  ring  shrinks  in  order  to  maintain  volume  unchanged.  After  sufficient  long  time,  the 
ring  collapses  into  a  ball.  During  this  time  evolution  process,  the  volume  of  the  cellular  construct 
is  preserved.  Figure  3  depicts  five  snapshots  of  the  simulation  at  t  =  0.1,2,16,30,38.  Compared 
with  the  fusion  into  a  smooth  ring  construct,  the  evolution  to  contract  into  a  solid  spherical  ball 
takes  much  longer  time  shown  in  Figure  3. 

4.3  Tissue  sheet  fabrication  via  single  layer  deposition  of  cellular  spheroids 

We  then  simulate  the  fabrication  of  a  tissue  sheet  bio-construct  via  a  single  layer  cellular 
spheroids  deposition.  We  lay  9  cellular  spheroids  initially  on  a  flat  surface  in  a  square  pattern. 
The  bio-construct  gradually  joins  into  a  network  construct  with  four  visible  holes.  The  entire  sheet 
keeps  contracting  to  close  the  holes  and  eventually  leads  to  a  solid  sheet.  The  time  scale  for  sealing 
the  holes  vis  fusion  is  longer  than  the  time  used  to  fusion  the  spheroids  in  contact  to  a  smooth 
networked  construct.  The  dissipation  of  the  cellular  material  and  the  surface  tension  are  the  major 
mechanisms  driving  the  fusion  process.  Figure  4  depicts  the  simulation  all  the  way  to  the  moment 
a  smooth,  solid  sheet  is  formed.  The  packing  of  the  initial  spheroids  is  not  optimal.  If  we  pack  the 
spheroids  tightly,  the  fusion  can  be  reduced  by  50%  or  more  as  we  will  shown  in  the  next  example. 

4.4  Lumenized  tube  generated  via  layer-by-layer  deposition  of  cellular  spheroids 

To  simulate  the  layer-by-layer  deposition  to  fabricate  the  vascular  vein  or  lumenized  tube,  we 
stack  rings  of  spheroids  one  on  top  of  the  other  up  to  three  layers  in  our  next  simulation.  The 
initially  rugged  tube-construct  evolves  into  a  smooth  tube  with  the  walls  thickened  over  time,  which 
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Figure  4:  A  tissue  sheet  is  fabricated  via  fusion  of  a  single  layer  deposited  cellular  spheroids.  The 
snapshots  are  taken  at  t  =  0, 0.5, 1, 2, 3. 


Figure  5:  The  tube-formation  in  layer-by-layer  deposition  of  cellular  spheroids.  Snapshots  are  taken 
at  t  =  0,0.1,0.24,0.25,0.5.  It  completely  fused  at  t  =  0.25. 

qualitatively  captures  the  morphological  evolution  of  the  fusion  process  observed  in  experiments 
(Figure  5). 

The  positioning  or  packing  of  the  spheroids  relative  to  each  other  seems  not  to  make  a  differ¬ 
ence  qualitatively.  Quantitatively,  however,  a  bioconstruct  with  an  initially  more  tightly  packed 
spheroids  can  accelerate  the  fusion/maturation  process  significantly.  Figure  6  depicts  fusion  of  a 
tightly  packed  spheroid  tube  construct,  in  which  the  tube  fuses  into  a  hole-less  tube  in  much  shorter 
time  than  the  loosely  packed  tube  construct  depicted  in  Figure  5. 

4.5  Bifurcating  vascular  junction  fabrication  via  layer-by-layer  deposition  of 
cellular  spheroids 

We  extend  the  layer-by-layer  deposition  strategy  to  fabricating  a  bifurcating  lumenized  tube.  If 
the  initial  deposition  of  the  construct  is  done  without  defect,  i.e.,  too  large  a  hole,  both  Y-shaped 
and  T-shaped  joint  can  be  formed  within  the  same  period  of  time  like  a  single  lumenized  tube. 


Figure  6:  The  tube-formation  in  layer-by-layer  deposition  of  cellular  spheroids.  Snapshots  are  taken 
at  t  =  0,0.008,0.009,0.2.  It  completely  fused  at  t  =  0.009. 
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Figure  7  depicts  a  Y-shaped  bifurcating  tube  fabricated  via  fusion  of  cellular  spheroids.  Figure  8 
portraits  a  T-shaped  bifurcating  tube  obtained  by  layer-by-layer  deposition  of  cellular  spheroids. 
If  the  initial  deposition  were  not  done  properly,  a  defective  bifurcating  tube  may  form  shown  in 
Figure  9,  where  permanent  holes  form.  Therefore,  the  initial  deposition  of  the  cellular  spheroids 
are  instrumental  to  ensure  the  quality  of  the  final  product.  This  numerical  tool  can  then  be  used 
to  study  the  correlation  between  the  initial  configuration  of  the  spheroid  arrangement  and  the 
possible  product  defect  due  to  the  unexpected  holes  on  the  wall.  Ultimately,  it  can  be  embedded 
in  a  computer-aided  design  software  for  the  bio-printing  fabrication. 


Figure  7:  Branching  vascular  vein  fabricated  via  bioprinting  technology.  Good  branching  vascular 
construct  made  of  layer-by-layer  deposition  and  strategic  remedy  for  the  initial  profile.  Snapshots 
are  taken  at  t  =  0, 0.2, 0.3, 0.5, 0.8. 


Figure  8:  T-shaped  branching  vascular  vein  fabricated  via  bioprinting  technology.  Perfect  branching 
vascular  vein  made  from  layer-by-layer  deposition.  Snapshots  are  taken  at  t  =  0, 0.2, 0.3, 0.5, 0.8. 


5  Conclusion 

We  have  developed  a  phase  field  model  for  studying  fusion  of  cellular  aggregates  by  treating  the 
cellular  aggregate  and  the  host  hydrogel  matrix  as  two  immiscible  fluids.  The  model  incorporates 
the  long-range  attractive  interaction  between  adjacent  cellular  aggregates,  philic  (attractive)  inter¬ 
action  in  close  contact  and  the  repulsive  interaction  due  to  immiscibility  of  distinctive  materials 
at  the  material’s  interface.  An  efficient,  high  order  spectral  method  is  employed  to  discretize  the 
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Figure  9:  Perforated  tube  made  via  layer-by-layer  deposition  of  cellular  spheroids  with  a  bad 
packing  method.  A  defective  branching  vascular  construct  forms,  in  which  holes  are  present  in  the 
branching  point.  Snapshots  are  taken  at  t  =  0, 0.2, 0.4, 0.6. 

governing  partial  differential  equation  in  the  model.  A  numerical  solver  is  developed  by  implement¬ 
ing  the  spectral  method  in  2-D.  A  3-D  numerical  solver  is  developed  based  on  the  Cahn-Hilliard 
equation  only  to  be  applied  to  investigate  the  fusion  process  in  which  the  cellular  spheroids  are 
placed  in  contact.  Both  solvers  are  then  applied  to  bio- fabrication  of  various  specific  geometric 
shaped  tissues  and  organs  via  layer-by-layer  deposition  of  spheroidal  cellular  clusters.  The  hydro- 
dynamic  effect  involving  momentum  transport  is  shown  to  be  important  in  the  initial  compaction 
process;  it  is  later  dominated  by  the  philic  interaction  among  the  cells  once  the  cellular  clusters  get 
in  contact  to  each  other.  The  time  for  fusion  of  the  cellular  clusters  depends  on  the  initial  packing 
of  the  cellular  spheroids.  Tight  packing  will  certainly  result  in  shorter  fusion  time.  The  bifurcating 
bio-constructs  also  depends  on  the  packing  at  the  turning  point.  A  computer  aided  design  tool 
is  therefore  necessary  to  virtually  simulate  the  entire  biofabrication  and  provides  guidance  to  the 
bioprinting  process  in  precision. 
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